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Abstract. We present two numerical simulations of an accretion flow from a rotating torus onto a compact object with and 
without a solid surface - representing a neutron star and a black hole - and investigate its influence on the process of jet 
formation. We report the emergence of an additional ejection component, launched by thermal pressure inside a boundary layer 
(BL) around the neutron star and examine its structure. Finally, we suggest improvements for future models. 
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1. Introduction 

Although jets are ubiquitous phenomena in many different as- 
trophysical objects, their formation is relatively unclear. We 
find jets in young stellar objects where they are driven by pro¬ 
tostars, in symbiotic stars (white dwarfs). X-ray binaries (neu¬ 
tron stars and stellar mass black holes) and active galactic nu¬ 
clei (supermassive black holes). The mass loss rate of all jets is 
found to be connected to the mass accretion rate of the underly¬ 
ing disk found in most objects (e.g. iLiy ioi 19971) . Therefore the 
necessary components seem to be well known and common to 
all objects. 

In jet formation models presented so far, the magnetic field 
seems to play a key role. The first analytical work study¬ 
ing magneto-centrifugal acceleration along mag netic field lines 
thread ing an accretion disk was done bv iRlandford & Pavnel 
dl982l) . They have shown braking of matter in azimuthal direc¬ 
tion inside the disk and its acceleration above the disk surface 
by the poloidal magnetic field components. Toroidal compo¬ 
nents of the magnetic field then collimate the flow . Numer ous 
semi-a nalytic models extended the work of lBlandford & Pavnel 
dl982h . which either were restricte d to self-similar solutions 
and their geometric limitations (e.g. IPudritz & Normanllll986t 

IVlahakis & Tsinganoslll998l Il999t iFerreira & C assell20 04l) or 

suggested non-self-similar solutions (e.g. Icamenzindlll99(l 

IPelletier & Pudritzil 992lfBreitmoser & CamenzincJ200()lf 

Another approach is to use time-dependant numerical 
MHD simulations to investigate the formation and collima- 
tion of jets. In most models, however, a polytropic equilib- 
rium accretion disk was regard e d as a b ounda r y condition (e.g. 
iKrasnonolskv. Li & Blandfordl ; 19991 12004 ;; I Anderson et al.l 
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120041 Goodson. Bohm & Wingle3ill999 i). The magnetic feed¬ 
back on the disk structure is therefore not calculated self- 
consistently. Only in recent years were the first simulations pre¬ 
sented including the accretion disk self-consistently in the cal- 
culations_of jet fo rmatio n (e. g.lCa sse & KennenJ2002H2004 

iKato. Mineshige & Shibatal2004l). 

IPringlel dl989t) proposed the idea of jet formation in the 
BL region and also in his_model _strong magnetic fields were 
driving the outflows. Torhell dl984h first assumed the lib¬ 
eration of energy in BL shocks to drive winds by thermal 
pressure. iLiviol (119991 ) pointed out that an additional source 
of energy beside the m agnetic one is needed to power jets. 
pTorbett & Gilden ( '19921) performed numerical simulations and 
found mass ejection only when they had not taken radiative 
cooling into account. However, their simulations were only 
one-dimensional for calculating the vertical structure of the 
BL. New examinations and modifications of this possibility of 
acc elerating plas ma close to the central object were done by 
ISoker & Regevi d2003l) involving SPLASHs (S Patio temp oral 
Localized Accretion SHocks) in the BL. Locally heated bub¬ 
bles expand, merge, and accelerate plasma to higher velocities 
than the local escape velocity. This scenario was introduced in 
analytic estimates. Now the numerical treatment needs to be 
improved using mu/fz-dimensional simulations - at first purely 
hydrodynamical ones, which are presented in this paper. 

In Sect. 13 we present our numerical simulations of an ac¬ 
cretion flow onto a neutron star with a solid surface and onto 
an accreting black hole without one, with which we investigate 
its effects on the jet formation process (Sect.0. Sections0and 
0 examine the structure of the accretion flow and of the addi¬ 
tional ejection component, while a discussion follows in Sect. 
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2. The numerical models 

In the following, we describe our computer code with equa¬ 
tions, the model geometry, and the parameters. 

2.1. The computer code 

With the code NIRVANA ( iZicder Um. Il999l) we solve the 
following set of differential equations of ideal non-relativistic 
magnetohydrodynamics 

dp 

—+ V (p v) - 0 
at 

—+ V(pv®v) = -Vp+-(VxB)xB-pVch 
at p 

de 

— —I- V (e v) = -p V v 
ot 

p = (y - l)e 
<9B 

— = Vx(vxB) (1) 

at 

( 2 ) 

with density p, velocity v, internal energy e , pressure p, mag¬ 
netic field B, gravitational potential <b, magnetic permeability 
p, and adiabatic constant y. 

A set of common boundary conditions, including inflow, 
outflow (open), mirror and anti-mirror and rotational symme¬ 
try, all with their usual meanings, has already been defined in 
NIRVANA. The code and its bo u ndary c onditions were tested in 
many simulations flZiegleri 199119991 . and references therein). 


2.2. Initial conditions 


Besides taking the more or less standard disk as initial condi¬ 
tion, another approach is to begin with a rotating torus inside 
the computational domain. One advantage of this setup is that 
all material is already inside the domain initially, so no matter 
source has to be implemented on the boundaries. In this case, 
we could use the standard boundary conditions of NIRVANA. 

Starting with the static hydrodynamics equations, from the 
momentum equation follows 



- VO + 



(3) 


where r is the cylindrical radius. With a poly tropic equation of 
state and the identity 

p y-1 1 

V(-) = - -Vp (4) 

p y p 

Eq. 0 can be integrated to 


7 P 
y- i p 


+ <h + 



dr' 


l(r') 2 

r' 3 


= W 0 . 


(5) 


Under the assumption of constant /, the integral can be solved 
to 


r 

y- 1 


p 1 r 

- + <&+-— = W 0 
p 2 r z 


(6) 


with which the density is then 


1 y-1 
k y 


Wo - O - 


1 T 

2 r 2 


i/(r-i) 


(7) 


The angular momentum of the torus is then dependent on its 
radial position Rq as 


KRo) = ^GMR q 


Ro 


Rq — 2 R s 


( 8 ) 


if a pseu d o-Newt onian gravitational 
dPaczvnski & Wiitall980l) 

GM 

< 1 > =- 

R-2R S 


potential 


(9) 


is chosen. R is the spherical radius, and all distances are now 
given in units of R g . The corresponding time scale to is then the 
inverse of the Keplerian period to ~ 10 4 s in our simulations. 
The density maximum of the torus was positioned to 8 R s . 
Inside the torus, the velocity components are then 


V <P 


KRo) 


v,. = Vfl = 0. 


( 10 ) 


Outside the torus, Keplerian rotation is assumed 


W = Ok r v r = Vg = 0 (11) 

with 

Ok =^GMR—?—\. (12) 

R - 2 K g r z 

To initialise the magnetic field and to assure that its diver¬ 
gence vanishes, we calculate the magnetic field components 
from the vector potential A defined as 


B = V x A. 


(13) 


The only considered component of A is set identical to the in¬ 
ternal energy 

A,p = e. (14) 


The initial magnetic field lines are then along isocontours of en¬ 
ergy and density, as we used a polytropic equation of state dur¬ 
ing initialisation. This results in the following magnetic field 
components 


Br = 


1 


R sin# 86 


— (sin0e) 


1 8 

B e = — —(Re) 
R8R 


B, = 0 (15) 


in spherical coordinates, or 


8 


1 8 


B r = -—e B z = - — (re) 
8z r 8r 


B d = 0 


(16) 


in cylindrical coordinates. Afterwards, these components were 
scaled to achieve an assumed plasma /i = p gas /p mag = 10 3 . 
Note that no global external magnetic field was implemented. 

We performed two simulations with different boundary 
conditions at the inner radial coordinate, one with open bound¬ 
aries (Run A) - describing a black hole - and one with anti¬ 
mirror conditions to model the solid surface of the central ob¬ 
ject (Run B) - a neutron star. The other boundary conditions 
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are open ones at the outer radial coordinate, rotation symmetry 
and anti-mirror symmetry at the inner and outer poloidal ( 9 -) 
coordinate, respectively, to simulate the rotational axis and the 
equatorial plane, and periodic conditions at both azimuthal (</>-) 
boundaries. The aim of these simulations was to investigate the 
influence of the solid surface. A third simulation (Run C) was 
performed similar to Run A, but in cylindrical coordinates to 
investigate the influence of different coordinate systems. Here 
the cylindrical radius starts at 2 R g ; i.e. a cylinder along the axis 
is cut out and the inner radial boundary is chosen to be open. 

3. The effect of a solid surface on the jet formation 
process 

3.1. The accretion and ejection components 

In Fig.0 the logarithm of density is plotted for the three runs 
after 45 fo. One can clearly see that activity is triggered much 
faster in Run B with the solid surface boundary. 

In Run A accretion sets in immed i ate ly be cause of the 
magneto-rotational instability iBalbus & Hawlevil 19981) inside 
the torus. After almost one revolution of the torus, it comes in 
contact with the central object, which leads to a re-distribution 
of matter in the now established disk. The inner part of the disk 
puffs up creating an expanding bubble. Along the rotation axis, 
a funnel - i.e. a region evacuated by centrifugal forces to den¬ 
sities that are three orders of magnitude below the surrounding 
- is created and its cross-section grows with time. 

In Run B the accretion process also starts immediately and 
the funnel is created. Additionally, a boundary layer with a ra¬ 
dial extent of one fifth of the inner radius forms on the surface 
of the central object. Due to its high pressure, a flare occurs. 
This flare bubble, now created by the boundary layer and not 
by the accretion disk, expands along the surface of the torus 
with high velocities and high density. The density contrast of 
the bubble is 77 = 10 2 , i.e. the bubble is overdense. Expansion 
of this bubble is powered by a continuous high pressure flow 
from the boundary layer. The expansion direction of the bubble 
is along a latitude of around 40-50°. If an external magnetic 
field were present, perhaps the flow would bend towards the 
axis. 

Only in Run C, a high velocity component emanates out of 
the spherical expansion inside the swept-out funnel along the 
axis after 64.5 fo- This funnel jet has a head velocity between 
c/3 and c and a density contrast of rj — 10~ 2 . As this is only 
seen in Run C and not in Run A, the open boundary might 
create this feature. 

In Fig. |3 the time evolution of the accretion rate in the 
equatorial plane at r — 4 R g is plotted. It is calculated as 
M - -4 up v r H with height H of the disk. The qualitative be¬ 
havior seems to be equal in all runs, but in Run C the accretion 
rate is always higher. It is possible that the open boundary near 
the axis causes a global radial inflow which is superimposed 
on the accretion rate visible in Runs A and B. The accretion 
rate in the Runs A and B seems to be equal until about 21 fo; 
i.e. no causal connection is established between the boundary 
layer and the accretion flow itself until then. After that point. 
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Fig. 1. Plots of the logarithm of density of Run A (top). Run B 
(middle) and Run C (bottom) after 45 fo; the coordinates are in 
units of R g 
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Fig. 2. Time evolution of accretion rate M — -Anpv r H in 
the equatorial plane at r — 4 R g . In Run C the accretion rate 
is always higher than in the other simulations, perhaps due to 
the boundary conditions. After 21 to, the flare coming from the 
boundary layer destabilizes the accretion flow and increases its 
rate in Run B. 


the flare coming from the boundary layer destabilizes the ac¬ 
cretion flow and increases its rate. 

3.2. Jet emission efficiency 

In Fig. 0 a poloidal slice at r — 10 R g of the mass accre¬ 
tion/outflow rate - again calculated as M — -4 npv r H - is 
plotted at different times for Run A and for Run B. One can 
clearly see the different accretion and ejection components in 
this simulation. 

The first peak at about 40° represents the flares created by 
the boundary layer in Run B and has no corresponding fea¬ 
ture in Run A. The second peak at about 70° is identical to an 
outflow along the surface of the torus, which is common in all 
runs. At larger values of 6, the accretion flow can be seen. At 
this distance of 10 R g the accretion rate is higher in Run A than 
in Run B, while at closer distance the behavior is reversed (Fig. 
|2}. The accretion rate in all peaks is highly time dependent in 
Run B, while it seems to reach an asymptotic value in Run A 
(Fig. a, which results from the flary conditions inside the BL. 

Using the mass accretion rate in the equatorial plane and 
the mass ejection rate of the two peaks, one can calculate an 
ejection efficiency of the system as the ratio of both rates. This 
efficiency is plotted in Fig. [3 After an initial phase of global 
ejection in the equatorial plane until about 14 to, the mass frac¬ 
tion outflowing along the torus surface (second peak) compared 
to that being accreted is almost constant between 25-30 % in 
Run A, while oscillating around a mean of about 50 % in Run 
B. In Run A, the efficiency of ejection in the first peak is only 
one percent and the ejection peak is not really present. In Run 
B, however, the ejection efficiency of the first peak is compa¬ 
rable to that of the second peak, i.e. also in the range between 
40-50 %. Therefore almost the whole accreted matter is ejected 
out of the central region. 




Fig. 3. Poloidal slice at r — 10 R g of the mass accretion/outflow 
rate M - -Anp v r H of Runs A (top) and B (bottom) at differ¬ 
ent times. The emergence of an additional outflow component 
is visible at about 40°. 

4. Structure of the accretion flow 

The next step is to investigate the structure of the accretion flow 
and to raise the question of influences of the emergence of the 
boundary layer. In Figs. l6l9( the main magnetohydrodynamical 
quantities of the flow are plotted. 

In Fig.^which shows the density distribution in the equato¬ 
rial plane, the boundary layer between 2 and 2.4 (1-1.2 radii 

of the central object) is, along with the torus, the most promi¬ 
nent feature in Run B. The density increases by more than an 
order of magnitude with respect to the accretion flow. In Run 
A, the open boundary representing the inner sink creates a den¬ 
sity decrease by a factor of 4-5, due to draining into the black 
hole. In contrast to Run B, in which the accretion flow is always 
rotating super-Keplerian, the angular momentum of the flow in 
Run A drops below the Keplerian angular momentum for dis¬ 
tances smaller than about 3 R g (Fig. 0- This is also an effect 
of the drag created by the black hole. The thermal pressure is 
enhanced in Run B inside the BL by six orders of magnitude 
compared to Run A (Fig. [3- The emergence of the boundary 
layer also creates a peak in the radial and azimuthal compo¬ 
nents of the magnetic field (Fig. [3 at its surface, which is per¬ 
haps caused by compression in strong shocks inside it. This 
highly magnetised region extends to a distance of 60 R g from 
the equatorial plane, leading to a positive total poloidal current 
instead of negative values near and inside the torus (Fig. mu. 
which then can also drive jets magnetically. 
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Fig. 4. Time evolution of the ejection peaks and the accretion 
disk in Run A (top) and Run B (bottom). While the accre¬ 
tion/outflow rates seem to reach an asymptotic value in Run 
A, they are highly variable in Run B due to the flary conditions 
inside the BL. 

5. Structure of the BL ejection component 

In Figs. 1 1 1 1 j~2l the density and temperature along a slice through 
the ejection component at 8 — 45° are plotted. The variability 
of the accretion and ejection rates leaves its mark in the form 
of a rich substructure in all variables, especially in density. One 
can distinguish knots with enhanced density, magnetic field 
components, and temperature. In the velocity components, the 
knots are modulated on a global deceleration of the flow. The 
jet is still highly transient and simulations with an extended 
computational domain have to show whether a steady outflow 
on larger scales can be established. 

6. Discussion 

We set up numerical simulations of a compact object with and 
without a solid surface accreting matter from a rotating torus. 
They show an additional ejection component that could be col¬ 
limated into a jet by a global magnetic field which was, how¬ 
ever, omitted in our simulations. Another possibility for achiev¬ 
ing collimation is an external pressure gradient, which was also 
not present in our simuations. 

Our results seem to support the SPLASH scenario of 
ISoker & Regevl d2003h . We reported the emergence of a ejec¬ 
tion component which is highly variable in agreement with 
their scenario. In our simulations a two-dimensional treatment 


Fig. 5. Efficiency of the ejection mechanism, calculated as the 
ratio of the outflow rate in the ejection peak to the accretion rate 
in the equatorial plane plotted for Runs A (top) and B (bottom) 



Fig. 6. Structure of the accretion flow after 14 to: density. The 
density increases in the BL by more than an order of magnitude 
with respect to the accretion flow, in Run A, the open boundary 
representing the inner sink arranges for a density decrease by a 
factor of 4-5. 

was used, the full description in three dimensions could reveal 
new effects or details which we have to look for in new simu¬ 
lations. 

Another point is the amount of physics in our models. We 
used the equations of ideal MHD neglecting any cooling ef¬ 
fects. Depending on the accretion rate, however, the bound¬ 
ary layer material can become optically thick, which has to be 
taken into account. The additional equations in the flux-limited 
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Fig. 7. Structure of the accretion flow after 14 to: ratio of angu¬ 
lar momentum to Keplerian angular momentum. In Run B the 
accretion flow is always rotating super-Keplerian and the angu¬ 
lar momentum of the flow in Run A drops below the Keplerian 
angular momentum for distances smaller than about 3 R g . 



Fig. 8 . Structure of the accretion flow after 14 to: pressure. The 
thermal pressure is enhanced in Run B inside the BL by six 
orders of magnitude with respect to Run A. 


diffusion ansatz can no longer be solved by NIRVANA, so that a 
new tool has to be used e.g. FLASH, and a new branch of simu¬ 
lations will be necessary. As cooling effects reduce the internal 
energy of the material and thereby also reduce its thermal pres¬ 
sure, it has to be shown, whether this new jet formation scenario 
is still reliable in two- or three-dimensional models. 

We have presented simulations with a set of model param¬ 
eters appropriate for an accreting neutron star and an accreting 
black hole, respectively. The accretion rate in our simulations 
is, however, far too large for XRBs, but would only be suitable 
for a gamma-ray burst. These are created by an overestimated 
density inside the rotating torus. Further simulations need to 
show whether this scenario still works at lower rates and will 
have to fix the val ue of a critical accretion rate, as stated in the 
analytic model hv ISoker & I .asotal ( ‘2004 4 

The equations of ideal MHD can theoretically be written in 
a non-dimensional form, if one uses a Newtonian gravitational 
potential instead of the pseudo-Newtonian one. Neglecting the 
latter, we could normalize all quantities to naturally arising 
combinations which depend on parameters of the central ob¬ 



Fig.9. Structure of the accretion flow after 14 fo: azimuthal 
magnetic field component. The emergence of the boundary 
layer also creates a peak in the magnetic field at its surface, 
which is perhaps caused by compression in strong shocks in¬ 
side it. 

ject and carry over our results to other jet sources. However, 
additional simulations representing other classes of jet emitting 
objects will follow. 
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Fig. 11. Structure of the ejection component at 6 = 45° after 54 
to: density. Knots with enhanced density are present. 
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Fig. 12. Structure of the ejection component at 6 — 45° after 54 
to: temperature 
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Fig. 10. Plots of the total poloidal current r B$ after 31 to for 
Run A (top). Run B (middle) and Run C (bottom) with loga¬ 
rithmic contour lines from / = 0 to / = 10 12 A 
















